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ABSTRACT 


The  fundamental  concepts  of  the  theory  of  invariant  imbedding 
as  applied  to  particle  transport  processes  are  reviewed  and  their 
applications  to  shielding  problems  are  discussed.  Two  new  approxi¬ 
mate  methods  for  predicting  the  transmission  of  particles  through 
slabs  of  material  are  derived  from  the  imbedding  equations.  One- 
velocity  computations  of  the  angular  distributions  of  particles 
reflected  and  transmitted  by  slabs  of  varying  thickness  and  com¬ 
position  are  presented  and  analyzed. 
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SUMMARY 


Invariant  imbedding  is  a  new  approach  to  the  derivation  of 
transport  equations  and  results  in  a  formulation  quite  different 
from  the  usual  Boltzmann  formulation.  It  yields  equations  for 
the  fluxes  transmitted  and  reflected  by  a  body  as  a  function  of 
the  dimensions  of  the  body.  For  this  reason  it  is  particularly 
well  suited  to  shielding  studies  in  which  shield  thickness  is 
considered  as  a  parameter.  This  report  summarizes  the  develop¬ 
ment  of  the  imbedding  theory  and  discusses  shielding  applications 
of  the  theory. 

The  imbedding  and  Boltzmann  approaches  are  compared  to 
illustrate  the  essential  differences  in  the  resulting  equations. 

It  is  noted  that,  in  general,  the  imbedding  method  yields 
initial-value  rather  than,  two-point  boundary-value  problems  as 
does  the  Boltzmann  method.  The  resulting  computational  advantages 
of  the  imbedding  formulation  are  discussed. 

The  imbedding  equations  for  the  transmitted  and  reflected 
energy  and  angular  distributions  in  a  slab  geometry  are  derived. 
Extension  of  the  method  to  other  geometries  and  multilayer  con¬ 
figurations  is  discussed.  Methods  for  applying  the  imbedding 
equations  to  the  prediction,  of  emergent  neutron,  gamma,  neutron- 
induced  gamma,  and  photoneutron  distributions  are  presented. 

It  is  shown  that  the  equations  for  these  processes  are  of  the 
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same  form  and^  further,  that  they  may  be  solved  simultaneously 
by  the  same  techniques  used  in  obtaining  primary  neutron  or  gamma^ 
ray  solutions.  The  use  of  imbedding  methods  in  shield  optimization 
studies  is  discussed  briefly.  Methods  for  obtaining  approximate 
solutions  to  thick— slab  penetration  problems  are  presented. 

A  FORTRAN  Procedure,  R61.  coded  for  the  IBM— 7090  to  solve  the 
imbedding  equations  for  a  slab,  is  described.  Numerical  results 
of  some  R61  calculations  of  one— veloci ty,  emergent  angular  distribu¬ 
tions  are  presented,  A  simplified  model  is  developed  and  applied 
in  the  analysis  of  these  data  to  facilitate  qualitative  comparisons 
of  angular  distributions  as  a  function  of  slab  thickness,  absorp¬ 
tion  probability,  and  angular  scattering  cross  section.  As  a 
result  of  the  numerical  studies  a  method  which  results  in  further 
Simplification  of  the  thick-slab  penetration  problem  is  developed 

It  is  concluded  that  the  method  of  invariant  imbedding  is 
potentially  a  powerful  tool  for  theoretical  studies  of  shielding 
problems  as  well  as  for  the  generation  of  parametric  data  for 
shield— design  studies. 
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I,  INTRODUCTION 


The  development  of  invariant  imbedding  began  with  the  work 
of  the  Russian  astrophysicist,  Ambarzumianj,  who  derived  an  equa¬ 
tion  for  the  reflection  of  light  from  an  infinite  half  space 
(Refo  1),  The  solution  was  obtained  by  Chandrasekhar  (Ref,,  2)„ 
Although  this  work  was  of  great  theoretical  importance  because 
it  constituted  a  new  approach  to  transport  theory,  it  was  largely 
ignored  by  workers  in  neutron  transport  theory  until  a  few  years 
age.  During  the  past  five  years  Bellman,  Kalaba  and  Wing  (Ref..  3) 
have  developed  a  generalization  of  Ambarzumian’s  idea  which 
includes  transport  equations  for  finite  media. 

Invariant  imbedding  is  not  a  method  for  obtaining  analytic 
or  numerical  solutions  to  the  Boltzmann  equation.  It  is  a  new 
approach  to  the  derivation  of  transport  equations  -  quite  different 
from  the  usual  Boltzmann  approach.  In,  fact,  it  is  more  general 
than  particle  transport  theory  itself.  Imbedding  techniques  are 
now  being  applied  to  problems  involving  wave  processes  (Ref,  3). 

For  our  purposes,  however,  it  may  be  considered  as  a  method  for 
obtaining  new  transport  equations  which  may  be  easier  to  solve 
than  the  Boltzmann  equation. 

In  the  following  sections  the  theory  and  some  possible 
applications  are  presented.  First,  the  imbedding  and  Boltzmann 
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equations  are  derived  for  a  simple  one^dimensional  case  and  the 
differences  between  the  two  approaches  are  discussed.  Next,  the 
imbedding  approach  is  applied  to  a  more  realistic  problem  that 
of  transport  through  a  slab.  Further  generalizations  are  dis¬ 
cussed,  The  last  two  sections  deal  with  the  application  of 
invariant  imbedding  to  shielding  problems  and  certain  approxima¬ 
tions  suggested  by  the  imbedding  approach. 

The  discussion  of  the  theory,  given  in  Section  II,  is  based 
primarily  on  notes  taken  during  a  course  in  transport  theory 
presented  by  G.  M,  Wing  at  the  University  of  California  at 
Los  Angeles,  The  subject  is  treated  in  much  more  detail  in  a 
forthcoming  book  by  Dr,  Wing, 
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II.  THEORY 


2 . 1  Comparison,  With  the  Boltzmann  Theory 

To  begin  the  presentation  of  the  theory,  both  the 
Boltzmann  and  imbedding  equations  are  derived  for  a  simple 
(admittedly  unrealistic)  problem.  The  purpose  at  this  point 
is  simply  to  illustrate  the  imbedding  approach  and  point  out  the 
differences  between  the  imbedding  and  Boltzmann  methods, 

2,1.1  Boltzmann  Method 

Consider  a  one— dimensional  rod  of  length  x  with  one 
particle./sec  incident  from  the  right  at  x.  No  particles  enter 
from  the  left.  Let  the  position  coordinate  in  (O, x)  be  z 
and  denote  the  left—  and  right-moving  fluxes  by  v(2)  and  u(2),. 
respectively..  Let  the  interaction  cross  section  be  a,  a 
constant;  and  let  F  and  B  be  the  probabilities  of  forward  and 
backward  scatter,  respectively. 

^  v(2) 

- ►u(z) 

6  ^  i+A  i  1>ne  particle/sec 

To  derive  the  Boltzmann  equations  for  u  and  v  consider  the 
fluxes  at  z  +  A . 
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u(z+A)  -  (l-oA)  u(4)  +  oAPu(z) 

+  aAB  v(z)  +  o(A)  ,  (1) 

where  the  terms  on  the  right  side  are,  to  order  A 9  the  flux 

due  to  particles  which  pass  from  z  to  z+A  having  no  interaction 

in  As  the  right-amoving  flux  that  scatters  forward  in  A,-  and 

the  left— moving  flux  scattered  backward  in  A.  The  term  0(Aj 

2 

contains  all  terms  of  order  A  or  higher.  One  might  think  of 
the  right  side  of  (1)  as  a  Taylor’s  series  expansion  about  z 
of  the  terms  describing  the  collision  processes  in  the  increment 
(z,  Z4-A),  Taking  the  limit  as  A— ►O,  we  have 


u(zfA)  -  u(z)  r  o(A>‘ 

lira  .  —  —  I  »  lim  -ou(  a)<»-oFu(z)-^oBv(z)  t - 

A-e  0  ^  A— ►  oL  ^  - 


du(  z) 
dz 


because 

•  (A) 

lira  -  ■  -  0 

A-^0  A 


<2> 
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Similarly, 


v(z)  -  (1—0 A)v(z+A)  +  oAFv(z+A) 

+  cABu(z)  +  o(A)  , 

v(z+A)  —  v(z)  o(A) 

- -  o(F-l)v(z+A)+oBu(z)  + - 

A  A 


dv(z)  _ 

- - =  a  j^(F-l)v(z)  +  Bu(z)J  ,  (3) 

dz 

The  boundary  conditions  are 


v(x)  -  1 

(4) 

u(o)  “  0  . 

(5) 

Equations  2  and  3  with  the  boundary  conditions  (Eqs.  4  and  5) 
constitute  the  Boltzmann  formulation  of  the  problem.  Note  that 
Equations  2  and  3  are  linear  and  that  Equations  4  and  5  are  condi¬ 
tions  at  the  boundaries  of  the  transport  region;  thus,  the 
arbitrary  constants  involved  in  the  solution  of  Equations  2  and 
3  are  determined  by  the  size  of  the  rod.  This  means  that  each 
different  size  rod  determines  a  new  problem.  If  one  wishes  the 
reflected  and  transmitted  fluxes  u(x)  and  v(0),  new  integration 
constants  must  be  determined  for  each  rod  length  before  these 
quantitl0s.«an  be  computed.  Of  course,  in  this  simple  example 
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the  Integration  constants  are  easily  determined  as  a  function  of 
X.  However,  in  more  complex  problems  the  determination  of  co:a~ 
stants  must  be  done  numerically,  thus  making  the  determination 
of  u(x)  and  v(0)  quite  complicated.  It  should  be  noted  that 
the  determination  of  the  integration  constants  of  a  system  of  N 
linear,  first-order  equations  requires  the  inversion  of  an 
N  X  N  matrix  for  each  boundary— value  problem.  (Systems  of 
coupled  equations  of  this  type  occur  in  nearly  all  solutions 
of  the  Boltzmann  equation;  e.g.,  Carlson  spherical  harmonics, 
and  double  P/, )  The  economics  of  matrix  inversion  on  a  computer 
is  such  that  these  equations  are  usually  integrated  numerically 
by  some  iterative  scheme.  This,  of  course,  does  not  avoid  the 
fact  that  each  new  configuration  presents  a  new  boundary-value 
problem  and,  thus,  a  new  calculation.  It  would,  therefore,  be 
highly  desirable  if  the  transport  problem  could  be  reformula tfd 
in  terms  of  the  dimensions  of  the  transport  medium.  The  > 

of  invariant  imbedding  provides  such  a  formulation  in  many 
cases. 

2,1.2  Invariant  Imbedding  Method 

Consider  once  again  the  rod  of  transporting  material  with 
one  particle/sec  incident  from  the  right.  The  right— moving  flux 
at  X,  u(x),  is  the  Reflected  flux  per  incident  particle/sec  which 
shall  be  renamed  R(x).  The  transmitted  flux  per  incident  particle, 
v(0),  shall  be  called  T(x).  (It  is  a  function  of  the  length,  x, 
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of  the  rod.)  The  problem  of  predicting  R(x)  and  T(x)  is  now 
"imbedded"  in  the  family  of  problems  of  predicting  R(x)  and 
T(x)  for  various  rod  lengths.  This  is  accomplished  by  consider¬ 
ing  the  changes  in  R(x)  and  T(x)  when  the  length  of  the  rod  is 
increased  from  x  to  x  -i-  A. 
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'■■■I 

x  x+A 


partlcle/sec 


Thus, 


R(x+A)  -  (l-oA)R(x)(l-aA)  +  oAB 
+  oAF  R(x)(l-aA) 

+  (l-oA)R(x)0^ 

+  (l-oA)R(x)oABR(x)(l-oA)+  o(A),  (6) 

where  R(x+A)  is  the  probability  of  reflection  by  a  rod  of  length 
x-i-A  and  R(x)  is  the  probability  for  length  x.  The  first  term 
on  the  right  is,  to  order  A ,  the  probability  that  a  particle 
will  suffer  no  collisions  in  going  from  x+A  to  x  times  the 
probability  of  reflection  from  the  "sub— rod"  of  length  x  times 
the  "no  collision"  probability  for  (x,^x+A).  The  second  term 
represents  the  back— scatter  probability  for  the  distance  from 
x+A  to  X.  The  third  term  represents  forward  scatter  in  (x,  x+A) 
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followed  by  reflection  from  the  sub-rod  with  no  subsequent  colli* 
sions.  The  fourth  term  represents  the  process  in  which  there  are 
no  collisions  in  (x,  x+  A)  followed  by  reflection  by  the  sub-rod 
followed  by  forward  scatter  in  (x,  x+A).  The  fifth  term 
represents  no  collision  in  (x,  x+A),  reflection  by  the  sub-rod, 
scatter  back  into  the  sub— rod  from  (x,  x+A),  another  reflection 
by  the  sub-rod  and  no  subsequent  collisions.  All  other  processes 

O 

are  of  order  A  or  higher  (i.e.,  they  involve  more  than  one  colli¬ 
sion  in  A)  and  are  therefore  included  in  o(A).  Rearrangement 
gives 


R(x+A)  -  R(x) 
A 


o 


B  +  FR(x)  -  2R(x)  +  FR(x) 


+  R(x)BR(x) 


o(A) 

A 


o  3 

where  0(A)  now  contains  the  A  and  A  terms  from  Equation  b. 
Letting  A — ►O,  we  have 


dR(x) 

dx 


B  +  2(F-l)R(x)  +  R(x)BR(x) 
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(7) 


a  first  order,  nonlinear  equation  for  R(x).  The  initial  condition 
is  R(0)  »  0,  which  simply  says  that  there  are  no  particles 
reflected  from  a  rod  of  zero  length. 
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similar  analysis  yields 


dT(x) 

- -  oT(x)(F-l  +  BR(x)), 

dx 


with  T(0)  =  1.  (8) 

At  first  glance  it  would  seem  that  the  problem  has  been  complicated, 
rather  than  simplified,  by  the  introduction  of  a  nonlinear  equation 
and,  it  must  be  admitted,  nonlinear  initial— value  equations  are 
typical  of  the  Imbedding  approach.  However,  one  must  remember 
that  virtually  all  practical  problems  in  transport  theory  require 
the  use  of  high-speed  computers  and  that  the  numerical  integration 
of  nonlinear  equations  is  no  more  difficult  than  the  numerical 
integration  of  linear  equations. 

The  additional  fact  that  Equations  7  and  8  represent  an 
initial— value  rather  than  a  two— point  boundary— value  problem 
makes  numerical  integration  easier.  Also,  because  these  are 
initial— value  problems,  it  is  not  necessary  to  carry  information 
concerning  all  mesh  points  in  the  fast  memory  of  the  machine;  only 
the  last  computed  values  are  needed  to  perform  the  next  step  in  the 
integration.  This  results  in  a  considerable  reduction  in  storage 
requirement,  which  is  usually  a  major  problem  in  transport  computa¬ 
tions.  In,  addition  to  these  computational  advantages,  the  imbedding 
formulation  provides  equations  with  the  material  dimension  as 
independent  variables.  Thus,  the  solutions  to  one  problem  are  the 
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transmission  and  reflection  values  as  a  function  of  the  physical 
dimensions  of  the  transporting  medium. 

2.2  The  Slab  Problem 

Consider  the  problem  of  energy-angular  dependent  transmis¬ 
sion  and  reflection  by  slabs.  To  simplify  the  algebra  the 
energy-angular  dependence  shall  be  represented  by  particle 
"states”  with  such  states  denoted  by  subscripts.  ThuSp  particles 
in  State  i  have  energies  which  lie  within  some  particular  energy 
group  and  are  traveling  in  a  particular  direction  (or  within  a 
group  of  directions).  Particles  in  the  same  energy  group  but 
with  different  directions  of  motion  are  considered  to  be  in  another 
state,  eay  i*.  The  choice  of  a  particular  mode  of  representation 
(such  as  discrete  angles,  angular  groups  or  Legendre  coef f icien it 
is  not  important  as  long  as  the  angular  dependence  can  be 
represented  by  discrete  numbers,  i.e.,  subscripts.  Consider 
slab  of  thickness  x  perturbed  by  the  addition  of  a  layer  of  thick¬ 
ness 
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Let 


o. 


3 


cross  section  for  interaction  of  State  j 
particles  with  the  inedium, 


=  probability  that  an  interaction  results  in 
the  transfer  of  particles  from  State  j  to  i 
through  a  forward  scatter, 

“  the  corresponding  probability  for  back— scatter, 
and 


Rj.  «  the  number  of  State  i  particles  reflected 
^  through  a  unit  area  on  the  surface  x  due 
to  one  State  j  particle/sec  Incident  per 
unit  surface  area. 


Note  that  the  probability  that  a  State  j  particle  will  suffer  a 
coll  ision  in  (x,x-f-A)  is  OjA/Wj ,  where  is  the  cosine  of  the 
angle  between  the  State  j  particle  direction  and  the  normal  to 
the  slab.  To  avoid  negative  values,  incident  angles  will  be 
measured  from  the  Inward  normal  and  reflection  angles  from 
the  outward  normal.  Proceeding  as  before,  it  is  seen  that 


(equation  continued) . 
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4- 


^  I 

k  JL 


AB,,  R  .U) 
/k  \l 


4  -  O  (  A  ,)  , 


where  the  sums  on  k  and  JL  are  over  all  particle  states, 
ing  and  taking  the  limit  as  before  one  obtains 


dx 


o 


ij 

J 


B 


R. . (x) 
ij 


3 


(X) 


/S 

F 

'ik 


,  I 

4 


Define  the  matrices 


F  - 


9 


where 


is  the  Kronecker  delta, 


B  - 


> 


Rearrange 
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and  R  -  (R^j ) . 

Then 

dR(x) 

- -  B  +  FR(x)  +  R(x)F  +  R(x)  BR(x).  (9) 

dx 

A  similar  development  gives 
dT(x) 

- -  T(x)  (F  +  BR(x)).  (10) 

dx 

The  initial  conditions  are 

R(0)  =  0  (null  matrix)  (11) 

T(0)  *  I  (identity  matrix)  (12) 

which  state  that  for  a  slab  of  zero  thickness  the  reflected  flux 
is  zero  and  the  transmitted  flux  is  the  same  as  the  incident 
flux;  i.e.,  Tij(O)  =  0  for  i  #j,  and  Tj^j^  =  1.  Equations  9  and 
10  are  a  system  of  first-order,  nonlinear  equations  with  initial 
values  given  by  Equations  11  and  12.  Numerical  solutions  may  be 
obtained  by  any  of  the  well-known  schemes  for  integrating  such 
systems  (Ref.  4).  A  FORTRAN  procedure,  R61,  has  been  programmed 
to  solve  Equations  9  and  10  by  the  Runge-Kutta  method  (see 
Section  IV).  Once  the  Matrices  R  and  T  have  been  determined, 


23 


transmitted  and  reflected  distributions  may  be  obtained  as  a  fu:nc«' 
tion  of  slab  thickness  for  any  incident  distributions.  The  pro« 
cedure  is  quite  simple.  Represent  the  incident  energy—angui a:r 
distribution,  by  intensities  in  discrete  states,  Sj .  Arrange 
these  numbers  in  a  column  matrix 


The  transmitted  distribution  for  a  slab  of  thickness  x  ,is  the 
column  matrix 


Similarly,  the  matrix  product  R(x)S  gives  the  reflected  distribu¬ 
tion,  Procedure  R61  computes  the  T  and  R  matrices  and  writes  these 
on  tape  for  any  number,  up  to  200,  of  designated  thicknesses.  It 
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is,  therefore,  possible  to  write  a  simple  matrix  multiplication 
procedure  which  would  read  T  and  R  data  from  the  R61  tape,  read 
incident  distribution  data  from  cards  or  tape,  and  compute  trans¬ 
mitted  and  reflected  energy— angular  distributions  as  a  function 
of  slab  thickness. 

Because  of  the  simplicity  of  the  slab  problem  it  was  the 
first  to  be  chosen  for  numerical  investigation.  Bellman, 

Kalaba  and  Prestrud  (Ref,  5)  have  performed  reflection  calcula¬ 
tions  for  a  one-velocity,  isotropic— scatter  model  for  slab 
thicknesses  up  to  20  mean  free  paths.  Their  computer  program  is 
similar  to  R61.  Bellman,  Kalaba  and  Prestrud  have  also  obtained 
time ->de pendent  results  for  a  model  with  no  energy  or  angular 
dependence  (Ref.  6),  During  checkout  of  R61,  the  author  ob¬ 
tained  transmission  and  reflection  data  which  are  presented  in 
Section  V 

2 . 3  Other  Geometries 

Imbedding  equations  for  cylinder  and  sphere  problems  have 
been  derived  by  Bellman,  Kalaba  and  Wing  (Ref.  3).  Numerical 
investigations  have  not  yet  been  attempted;  Bellman  and  Kalaba 
are  planning  to  run  some  sphere  problems  in  the  near  future. 

Wing  has  developed  a  generalized  approach  to  the  derivation  of 
equations  for  other  geometries  (Ref.  7)  and  has  investigated 
methods  for  transforming  from  the  Boltzmann  formulation  to  the 
imbedding  equations. 
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Multilayered  configurations  offer  no  special  difficulties.. 
Position— dependent  cross  sections  may  be  incorporated  in  the 
derivation  without  changing  the  usual  procedure.  The  slab  equa¬ 
tions  with  variable  cross  section  are  the  same  as  Equations  9 
and  10  except  that  the  F  and  B  matrices  are  functions  of  x. 

A  simpler  procedure  for  handling  multilayer  slabs  involves 
changing  the  initial  conditions.  Consider  a  rod  of  two  materials. 
Let  the  region  from  0  to  y  be  of  one  material  and  the  region 
from  y  to  the  right  be  another  material.  It  may  be  shown  that 

(5  y  X4-y 

the  reflection  solution  for  the  composite  rod,  x  -v  ni-iy  be 
obtained  by  solving  for  the  reflection  from  the  homogeneous  rod 
(0,  y)  and  then  proceeding  with  the  solution  for  the  second 
homogeneous  rod  (y,  x-j-y)  subject  to  the  new  initial  condition 
R(x=0)  =  R(y),  A  similar  modification  of  the  initial  condition 
may  be  used  in  solving  for  transmission  functions. 

This  procedure  may  be  adapted  to  any  multilayer  configura¬ 
tion. 
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III„  SHIELDING  APPLICATIONS 


3,1  The  Combined  Neutron  and  Ganima~Ray  Problem 

The  application  of  imbedding  techniques  to  neutron  and  gamma- 
ray  penetration  and  reflection  problems  follows  directly  from  the 
development  outlined  in  the  preceding  sections.  Application  to 
neutron»induced  gamma  calculations  is,  perhaps,  not  so  obvious  in 
view  of  the  fact  that  the  imbedding  method  does  not  provide  pre— 
dicticns  of  fluxes  inside  the  material.  It  shall  be  shown  in  this 
section  how  the  method  can  be  used  not  only  to  produce  secondary- 
gamma  predictions  without  explicit  knowledge  of  the  internal  neu¬ 
tron  flux,  but  also  to  calculate  transmission  and  reflection  of 
neutrons,  secondary  gammas  and  primary  gammas  simultaneously. 

In  Section  2,2  energy— angular  distributions  were  represented 
by  dividing  the  energy-angular  intervals  into  discrete  regions 
or  "states.”  Cross  sections  and  transfer  probabilities  were 
defined  for  each  state.  As  far  as  the  mathematical  development 
which  followed  is  concerned,  there  is  no  reason  to  assume  that 
particles  in  different  states  are  even  the  same  kind  of  particles. 
The  equations  are  still  valid  if  the  particles  in  the  first  few 
states  are  neutrons  and  those  in  the  remaining  states  are  gamma 
rays.  To  illustrate,  consider  a  simple  two— state  problem  in 
which  particles  in  State  1  are  neutrons  and  those  in  State  2  are 
gamma  rays.  Then  and  kre,  respectively,  neutron  and 
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gamma-ray  collision  cross  sections,  ®11!>  ^22  ®22 

neutron  and  gamma— ray  forward- and  back-scatter  probabilities, 
r2j^  is  the  probability  that  a  neutron  collision  results  in  a 
gamma— ray  traveling  in  the  precollision  neutron  direction,  Bgj 
has  a  similar  definition^  but  implies  a  reversal  of  direction, 

FjL2  and  gamma-to-neutron  conversion  and  will  therefore 

be  zero  unless  photoneutron  reactions  are  considered.  The  exten¬ 
sion  to  more  than  two  states  (to  allow  energy-angle  considerations 
to  be  incorporated)  is  straightforward  and  results  in  partitioning 
of  the  F,  B,  R  and  T  matrices  as  illustrated  below. 


The  sub-matrix  (n - m)  contains  the  probabilities  of 

transfer  from  one  neutron  state  to  another^  (n — contains 
neutron— to— gamma  probabilities,  etc.  Thus,  the  transmission  matrix 
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coTutains  the  transmission  matrices  for  neutrons,  primary  gamma 
rays,  secondary  gamma  rays  and  photoneutrons.  The  source  dis¬ 
tribution,  to  which  the  transmission  operator  is  applied  is 


where  is  the  column  matrix  representing  the  incident  neutrons 
and  S^.  represents  incident  gamma  rays.  The  result  of  the  opera¬ 
tion  “  TS  is  a  column  matrix  containing  the  transmitted  neutrons 
and  gamma  rays,  including  secondary  gamma  rays  and  photoneutrons. 
It  should  be  noted  that  the  equations  for  the  combined  neutron 
gamma=iay  problem  are  of  the  same  form  as  those  for  single- 
particle  transport.  Thus,  the  existing  code,  R61,  can  be 
applied  to  the  combined  problem. 

3.2  Shield  Optimization 

The  imbedding  method  may  also  be  helpful  in  shield  optimiza¬ 
tion  studies.  The  gradient-nonlinear  programming  method  (Ref.  8) 
now  used  to  arrive  at  an  optimum  shield  configuration  requires  a 
knowledge  of  the  derivative  of  the  dose  rate  (ideally  the  total, 
neutron  plus  gamma,  dose  rate)  with  respect  to  shield  weight  or 
thickness.  This  is  presently  accomplished  by  solving  shield 
transmission  problems  for  several  shield  thicknesses  (neutron. 
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gamma  and  secondary-gamma  predictions  require  separate  treatment), 
and  then  curve— fitting  and  differentiating  the  results.  An 
imbedding  formulation  of  the  same  problem  would  yield  a  differ¬ 
ential  equation  for  the  transmission  function  for  each  segment 
of  the  shield.  These  segments  might  be  considered  sections  of  a 
sphere,  in  which  case  the  imbedding  equations  already  exist  if  one 
is  willing  to  consider  each  segment  independent  of  all  others.  (This 
assumption  is  made  in  current  shield— penetration  procedures.)  The 
imbedding  equations  give  the  flux,  or  dose— rate,  derivatives  with 
respect  to  thickness  explicitly  in  terms  of  the  transmission  and 
reflection  functions.  Thus,  the  entire  optimization  procedure, 
with  the  exception  of  the  air-transfer  operation,  may  be  reformulated 
in  terms  of  transmission  and  reflection  matrices.  As  noted  in  the 
preceding  section,  these  matrices  may  include  all  processes, 
including  secondary-gamma  production,  and  may  be  obtained  for  a 
range  of  shield  thicknesses  in  a  single  problem.  Also,  note 
that  multilayer  configurations  are  easily  handled  as  long  as  thick¬ 
ness  variations  are  performed  with  the  outer  layer  only. 

3. 3  Approximations  for  a  Slab 

Note  that  as  the  thickness  of  a  slab  increases  the  reflec¬ 
tion  matrix  approaches  a  constant  matrix,  namely,  the  reflection 
matrix  for  an  infinite  half  space,  R^.  Thus,  for  x  large  enough, 
the  transmission  matrix  satisfies 
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dT(x) 


- -  T(x)(F  +  BR  ). 

dx  A 


The  solution  is 

(F+BR^)x 

T(x)  ••=  e  , 

where  the  well-known  matrix  exponential  is  defined  by 

/A  \  2  3 

(A)x  X  X 

e  =  I  +  Ax  4-  A'*  - +  A"^ - +  .  .  .  . 

2!  3! 

The  constant  matrix  is  the  solution  to  AJibarzumian ' s  problem 
and  may  be  obtained  by  Chandrasekhar’s  method  or  by  numerically 
integrating  the  imbedding  equation  until  sufficiently  small 
changes  in  R<x)  are  noted  at  neighboring  points.  The  matrix 
exponential  may  be  computed  either  from  the  definition,  from 
Sylvester’s  theorem,  or  by  diagonal! zation  of  the  matrix  (F  +  BR^). 

The  temptation  to  approximate  T(x)  by  a  matrix  exponential 
leads  to  another  approximation.  Note  that  the  term  BR(x}  represents 
a  double  reflection  because  B  is  the  back— scatter  matrix  and  R  is 
the  reflection  matrix.  If  scattering  is  predominantly  forward, 
one  may  wish  to  Ignore  multiple  reflections  of  the  type  illustrated 
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bp low  and 


represented  by  the  term  BR,  It  is  not  necessary  to 
assume  that  either  B  or  R  is  negligible,  onxy  viiat  the  product  BR 
is  negligible  compared  to  the  forward-scatter  matrix  F.  Then  the 
transmission  is  approximated  by 

Fx 

T(x)  =■  e 

and  the  reflection  by  the  solution  of 

dR(x) 

■■  B  FR  RF  • 

dx 

In  the  reflection  equation  the  term  RBR  has  been  neglected. 
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IV.  NUMERICAL  STUDIES 


A  series  of  one-veloeity  slab  penetration  and  reflection  com¬ 
putations  were  performed  that  will  illustrate  the  application  of 
the  imbedding  method  to  shielding  studies,  A  brief  discussion 
of  the  FORTRAN  Procedure  R61  is  given  below,  followed  by  a  detailed 
discussion  of  the  computations.  Some  of  the  computed  data  are 
presented  aiong  wfith,  an  analysis  based  upon  a  simplified  model  for 
describing  angular  distributions. 

4.1  FORTRAN  Procedure  R6X 

A  FORTRAN  Procedure,  R61,  was  coded  for  the  IBM-704  to  investi¬ 
gate  the  transmission  and  reflection  matrix  equations.  The  code 
is  based  on  the  Runge-Kutta  method  (Ref.  4)  and  includes  a  check 
of  the  results  after  each  step  in  the  integration.  A  description 
of  some  of  the  features  of  the  procedure  follows. 

In  general,  the  numerical  integration  of  an  initial— value 
problem  begins  with  a  knowledge  of  the  initial  value  of  the 
dependent  variable,  f<x),  and  its  derivatives, and  proceeds,  through 
a  series  of  finite  steps  in  the  independent  variable,  to  predict 
the  value  of  the  dependent  variable  at  the  end  point  of  each  step. 
The  Runge-Kutta  method  provides  a  formula  for  performing  this 
stepping  operation  in  terms  of  the  function  and  its  first  deriva¬ 
tive  at  the  last  mesh  point  considered.  That  is,  the  procedure 
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for  computing  the  number  f(x+A)  in  terms  of  f(x)  and 


X 


is  given  by  the  Runge~Kutta  formula. 

The  accuracy  of  the  computation  depends,  of  coarse,  on  the 
size  of  the  step  A  and  the  behavior  of  the  function  f(x)  in  the 
interval  (x,x+A).  An  estimate  of  the  error  involved  in  the 
computation  may  be  obtained  by  first  computing  f(x+A)  using 
the  full  step  A,  then  recomputing  f(x+  A)  by  taking  two  succes¬ 
sive  steps  of  width  A/2,  The  error  in  the  second  (more  accurate) 
estimate  of  f(x+  A)  is  approximately  one— fifteen th  the  difference 
between  the  first  and  second  computations  of  f (xf  A)  (Ref.  4). 

This  procedure  is  used  in  R61  to  estimate  the  error.  If  the 
error  in  the  prediction  of  f(x-4-A)  is  larger  than  some  pre— 
assigned  fraction  of  f (x-f  A)  the  entire  process  is  repeated  with 
A  replaced  by  A./2.  This  iterative  process  is  repeated  until 
either  the  error  is  small  enough  or  until  a  preassigned  maximum 
number  of  iterations  have  been  attempted.;  in  the  latter  case, 
the  problem  is  terminated.  The  computer  time  required  for  this 
checking  procedure  is  significant  and  hence  undesirable  in  a 
production  code.  It  was  incorporated  in  R61  to  provide  a  means 
of  studying  the  effect  of  step  size.  Some  of  the  results  of 
these  studies  are  given  in  the  next  section. 
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Because  the  Runge-Kutta  method  requires  the  function  and  its 
derivative  at  the  last  mesh  point  only,  it  was  possible  to  save 
space  in  the  fast  memory  of  the  machine  by  printing,  writing  on 
tape,  or  erasing  data  not  pertaining  to  the  last  mesh  point.  This 
reduced  storage  requirement  allows  consideration  of  more  energy 
or  angular  groups  than  would  ordinarily  be  possible. 

The  forward- and  back-scatter  matrices,  [^ijl  L^ij] 

(Sec.  2.2),  are  computed  in  the  program  from  the  formulas 


.  <5'  .. 

1  ’  J '  gg 
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L  2i+  1 
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2J^  ^  1 


i  =  MW  (g-l>  -  i^. 
j  »  MW  (g'-l)  4  j', 

where 

Og  is  the  total  cross  section  for  energy  group  g, 
u/j^t  Is  the  cosine  of  the  angle, 

•'ij  Is  the  K.ronecker  delta. 


Qj^’  is  the  weight  factor  associated  with  angle  i' 
arising  from  the  representation  of  an  integral 
over  angle  by  a  finite  sum, 
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is  the  Jith.  Legendre  coefficient  of  the  angular 
group  transfer  cross  section  for  scattering  from 
energy  group  g'  to  g  (Ref.  9)s 

is  the  ^th  Legendre  polynomialp  and 

is  the  number  of  discrete  ordinates  considered  in 
the  angle  quadrature. 

The  weights,  occur  when  the  imbedding  equations  are  derived 

with  continuous  angle  dependence  and  then  reduced  to  the  discrete 
form  indicated  in  Section  2.2  by  a  suitable  quadrature  formula. 

In  this  case,  the  first  term  on  the  light  side  of  Equal  ton.  9 
(Sec,  2.2)  is  not  the  B  matrix  defined  above,  but  another  matrix, 
B*,  defined  by 


■  * 

1  L  2^  +  1  .  . 

“ij 
..  • 

c  /:o  2  /  j 

The  derivation  of  these  relations  from  the  definition 


4' 


MW 


O  ,  4,  j 
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2t^  1 


4rr 


gg 


Pj(w) 


is  straightforward.  Here  Oggf(w)  is  the  cross  section  for  scat¬ 
tering  from  group  g’  to  g  via  the  scattering  angle  cos”^u . 

The  coefficients,  ,,  which  are  input  for  R61,  are  computed  in 

gg 

another  program,  C54  (Ref.  9). 
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The  program  output  consists  of  the  transmission  and  reflec¬ 
tion  matrices,  printed  and  written  on  tape,  at  preselected  slab 
thicknesses.  The  printed  output  is  a  square  array  of  numbers, 
each  column  ( j )  of  which  represents  a  particular  initial  state 
and  each  row  (i)  of  which  corresponds  to  a  particular  final  state. 
Thus,  the  number  located  in  the  array  according  to  conven- 

tional  matrix  notation,  is  the  probability  of  reflection  or 
transmission  for  initial  State  j  and  final  State  i.  The  States 
i  and  j  are  related  to  the  initial  and  final  energies  and  angles 
by  the  relations  given  earlier. 

4,2  Angular  Distribution  Computations 

The  cal  eolations  presented  in  Section  4.3  are  based  on  a 
five-point  Legendre-Gauss  approximation  to  the  angle  integration. 
Thus,  the  are  the  Legendre-Gauss  weights  and  the  are 
the  corresponding:  arguments  for  the  interval  (0,  1),  Because 
the  problems  involve  no  energy  degradation,  the  coefficients 
Aggt  reduce  to  the  Legendre  coefficients  in  the  expansion  of  the 
angular  cross  section.  For  simplicity  the  total  cross  section 
was  set  equal  to  unity  which,  in.  effect,  is  the  same  as  measuring 
all  distances  in.  mean  free  paths.  The  spatial  mesh  points  were 
so  chosen  that  the  interval  widths  become  larger  as  the  slab 
thickness  increases.  This  was  intended  to  take  advantage  of 
the  fact  that  the  most  rapid  changes  in  both  the  transmission  and 
reflection  functions  occur  for  small  thicknesses.  As  mentioned 
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earlier  (Sec,  4,1),  R61  will  reduce  the  width  of  an  interval 
if  the  error  involved  in  that  step  is  larger  than  some  pre¬ 
assigned  fraction  of  the  function.  For  the  problems  presented 
here  a  maximum  error  of  1%  per  step  was  allowed.  In  one  of 
the  problems  the  intervals  were  chosen  as  large  as  one  mean 
free  path  near  the  end  of  the  integration,  to  allow  the  program 
to  determine,  within  a  factor  of  2,  the  optimum  interval  for 
that  thickness  (near  10  mean  free  paths).  Three  interval 
divisions  occurred  which  left  the  interval  width  at  1/8  mean 
free  path.  Because  the  reflection  function  was  changing  very 
little  at  that  thickness,  the  interval  requirement  is  believed 
to  have  been  determined  by  the  transmission  function  only.  This 
indicates  that  if  reflection  data  are  all  that  are  desired, 
the  omission  of  the  transmission  calculation  could  reduce  the 
machine  time  by  more  than  a  factor  of  2  through  the  use  of 
larger  mesh-point  intervals  for  large  thicknesses. 

To  begin  the  integration,  step  widths  of  about  0,010  mean 
free  path  were  necessary  for  thicknesses  up  to  about  0.1  mean 
free  path.  Beyond  that  point  it  was  possible  to  Increase  the 
intervals  gradually  up  to  about  0,10  mean  free  path.  Apparently, 
when  both  transmission  and  reflection  calculations  are  to  be 
performed,  it  would  be  more  economical  to  use  a  constant  interval 
width  after  the  first  few  steps  and  thus  avoid  the  checking 
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procedure  after  each  step.  The  proper  constant  interval  could  be 
determined  by  using  the  interval  division  process  out  to  some 
pointj,  say  1  mean  free  path,  and  holding  the  Interval  constant 
from  that  point  on.  Elimination  of  the  check  after  each  step 
should  cut  computing  time  requirements  to  2/3  that  of  R61. 

The  choice  of  the  Runge—Kutta  method,  for  the  integration 
was  rather  arbitrary.  It  may  be  possible  to  reduce  the  computer 
time  by  using  other  standard  techniques  or  by  developing  special 
quadrature  formulas.  The  asymptotic  behavior  of  the  reflection 
and  transmiss.ion  matrices  might  be  used  as  a  guide  in  developing 
such  procedures. 

Computer  times  for  the  three  R61  problems  presented  in  the 
next  section  averaged  10  minutes  each.  The  number  of  mesh 
points  required  per  problem  varied  from  about  200  to  250. 

4.3  Analysis  of  Data 

The  numerical  results  presented  in  this  section  were  ob¬ 
tained  during  checkout  of  R61.  The  problems  were  chosen  to 
provide  data  for  an  analysis  of  the  angular  distribution  of 
particles  reflected  and  transmitted  by  slabs.  The  analysis 
was  restricted  to  one— velocity  problems  and  simple  angle- 
dependent  cross  sections.  While  the  generation  of  data  for 
practical  problems  is  possible  with  R61,  it  is  expected  that  a 
full-scale  attack  at  this  time  would  yield  only  a  multitude 
of  numbers,  overwhelmingly  complicated  by  energy-angle 
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relationships ;  uncertainties  in  slowing-down  approximations,  and 
the  usual  numerical  difficulties  encountered  in  manipulating 
large  systems  of  equations.  It  was  expected  that  the  study  of 
simple  problems  would  lead  to  a  better  understanding  of  the 
computational  and  physical  processes  involved  and,  thus,  to  more 
efficient  methods  for  realistic  problems. 

The  analysis  which  follows  will  therefore  be  restricted  to 
qualitative  discussions  of  the  effects  of  angular  cross  section, 
absorption,  and  slab  thickness  on  emergent  angular  distributions. 
These  discussions  will  be  based  on  certain  considerations 
of  the  collision  density  within  the  slab. 

4.3,1  Method  of  Analysis 

The  Integral  form  of  the  one— velocity  Boltzmann  transport 
equation  for  a  slab  is 


f (z, w ) 


S-0 


f(z',w')  o(A’ — ►A)dA’  d3  , 


where  f(z,w)  is  the  angular  flux  at  z  in  the  direction  A. and  S 
is  the  slant  distance  illustrated  on  the  following  page. 
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z— z* 

The  substitution  S  - -  gives 


i(z,w) 


f  (z*  ,g’  I^^)d2l*dz’ 


for  w]>0.  The  equation  for  the  number  of  particles/sec  passing 
through  a  unit  area  parallel  to  the  slab  faces  is  obtained  by 


multiplying  by  w 


uf  (z,iJ) 


z  z— z’ 

e  ^  (n 


f  (z'  )aCA'— •in.)d/\.'dz’  . 


n. 
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The  mean  value  theorem  gives 


uf  (z,u) 


(j 


f  (z*  ,(j*  )o(/t» — dz*  , 


0  Jv 


where  £(0  5^<z)  is  some  average  attenuation  distance  determined 
by  the  collision  density,  I  f  (z”  ,(j  '  )a(^ - feTv' )dA.«  ,  on  the 

■k 

conical  surface  defined  by  u  (see  sketch). 


0  Z'  z 


If  the  flux  within  the  slab  has  a  large  maximum  near  z*  »  0, 
the  average  attenuation  length  will  be  close  to  z.  On  the  other 
hand,  if  the  flux  is  rather  flat  as  a  function  of  z*  the  average 
length  will  be  less  than  z.  In  the  first  case,  the  transmitted 
angular  distribution  at  z  will  tend  to  be  more  strongly  peaked 
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toward  the  slab  normal  than  in  the  second  case.  This  can  be  seen 

z-$ 

by  examining  the  factor  e^^T  as  a  function  of  u  for  5-0 


and  ^  ]>  0.  Of  course,  the  angular  dependence  of  the  collision 

density  also  affects  the  resulting  distribution,  but  in  a  more 
complicated  manner.  This  effect  will  be  considered  briefly  in 
the  numerical  study  to  be  presented  later. 

Reflected  angular  distributions  may  be  studied  by  consider¬ 
ing  the  equation  for  u  0  (see  sketch): 

X 


L|f(z,  (j) 


f  (z* ,  (J  )o(3\-*— ►i^)d JT* 


dz' 


z  "A- 


where 

z  ^  4  <  X. 


43 


In  this  case  a  rapidly  decaying  spatial  distribution  implies 
^  ^  z  which,  in  turn,  implies  that  the  attenuation  term  is  of 
little  significance  except  for  small  u  .  It  would  then  be 
expected  that  the  angular  collision  density  would  be  the  dominant 
factor  in  determining  the  shape  of  the  reflected  distribution.  If, 
however,  the  flux  distribution  did  not  depend  so  strongly  on  z”, 
the  increased  importance  of  the  attenuation  term  would  be  expected 
to  cause  some  peaking  of  the  reflected  distribution  toward  the 
slab  normal. 

Although  the  preceding  arguments  were  applied  to  internal 
angular  fluxes,  they  are  also  applicable  to  surface  distributions 
and  may,  therefore,  be  used  to  analyze  reflected  and  transmitted 
distributions  as  a  function  of  slab  thickness.  If  the  internal 
flux  distribution  is  strongly  peaked  at  one  surface  of  the  slab, 
an  increase  in  slab  thicknesses  should  have  little  effect  on  the 
shape  of  either  the  reflected  or  transmitted  angular  distribution 
except,  perhaps,  for  some  increased  peaking  in  the  transmitted 
distribution.  Of  course,  a  large  variation  in  the  internal  flux 
from  one  surface  to  the  other  implies  a  rather  thick  slab  to 
begin  with.  If  there  is  little  variation  in  the  flux  throughout 
the  slab,  as  would  be  the  case  in  a  thin  slab,  an  increase  in 
thickness  would  be  expected  to  produce  more  peaking  in  both  the 
reflected  and  transmitted  distributions. 
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To  determine  the  significance  of  the  various  effects  Just 
described,  it  was  necessary  to  perform  calculations  of  reflected 
and  transmitted  fluxes  as  a  function  of  slab  thickness,  absorp¬ 
tion  and  angular  scattering  cross  section.  Another  parameter  of 
some  interest  is  the  incident— particle  direction  as  it  influences 
the  internal  flux  distribution  within  the  first  few  mean  free 
paths  from  the  entrance  boundary.  Information  concerning  the 
effects  of  slab  thickness  and  incident  direction  is  obtained 
from  a  single  problem  when  the  imbedding  method  is  used  because 
they  are  the  independent  variables  in  the  imbedding  formulation. 

To  study  the  effects  of  absorption  and  angular  cross  section  it 
was  necessary  to  run  additional  problems.  Three  cases  were  con¬ 
sidered;  1)  isotropic  scatter  with  the  absorption  cross  section 
equal  to  1%  of  the  total  cross  section,*  2)  isotropic  scatter 
with  the  absorption  cross  section  equal  to  one-half  the  total 
cross  section;  and  3)  anisotropic  scatter  with  no  absorption. 

The  anisotropic  cross  section  used  was  that  for  hydrogen: 
o(  cj)  ■■=  ij /it  for  w  >0,  c(w^  «  0  for  u<C0,  where  w  is  the 
angle  of  scatter.  The  data  presented  in  this  section  are  in  terms 
of  particles/sec  passing  through  a  unit  area  on  the  slab  surface 
due  to  one  particle/sec  entering  a  unit  area  on  the  surface. 

All  distances  are  measured  in  mean  free  paths.  Note  that  the 
source  strength  is  given  in  terms  of  particles  passing  through  a 
unit  area  on  the  surface;  the  data  should  therefore  be  multiplied 
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2  TT  X  ANGULAR  FLUX  (parficles/sec-unit  area-steradian) 


NPC  14,407 


0  2  4  6  8  10  12  14 

SLAB  THICKNESS  (mean  free  paths) 


FIGURE  1.  REFLECTED  ANGULAR  FLUX  AS  A  FUNCTION  OF  SLAB 
THICKNESS!  ABSORPTION  PROBABILITY  =  0.5 
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2  TT  X  ANGULAR  FLUX  (parficles/sec- unit  area-steradian) 


0  20  40  60  80  100  120  140 

REFLECTION  ANGLE  (deg) 

FIGURE  2A.  REFLECTED  ANGULAR  FLUX  AS  A  FUNCTION  OF  REFLECTION 
ANGlEi  ADSORPTION  PROBABILITY  =  0.5 
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2  TT  X  ANGULAR  FLUX  (particles/sec-unit  area-steradian) 


NPC  14,409 


0  20  40  60  80  100  120  140 

REFLECTION  ANGLE  (deg) 

FIGURE  2B.  REFLECTED  ANGULAR  FLUX  AS  A  FUNCTION  OF  REFLECTION 
ANGLE:  ABSORPTION  PROBABILITY  =  0.5 
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2  n  X  ANGULAR  FLUX  (par+icles/sec-unit  area-steradian) 


NPC  14,410 


0  2  4  6  8  10  12  14 

SLAB  THICKNESS  (mean  free  paths) 

FIGURE  3.  REFLECTED  ANGULAR  FLUX  AS  A  FUNCTION  OF  SLAB 
THICKNESS:  ABSORPTION  PROBABILITY  =  0.01 
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by  the  cosine  of  the  incident  angle  to  convert  to  a  source  of  one 
particle /sec  per  unit  area  normal  to  the  particle  direction, 

4,3,2  Reflection  Calculations 

Figures  1,  2A,  and  2B  show  the  reflected  angular  flux  for 
isotropic  scattering  with  the  absorption  cross  section  equal  to 
one<"half  the  total  cross  section.  The  first  figure  shows  that 
for  slab  thicknesses  greater  than  about  3  mean  free  paths  the 
reflection  is  essentially  constant.  For  other  incident  angles 
(not  shown)  the  behavior  is  similar.  Figures  2A  and  2B  show 
that  the  reflected  angular  distributions  for  1,1~  and  10,l~mean>^ 
free^-path  slabs  are  about  the  same.  These  distributions  are 
rather  flat,  especially  as  the  angle  of  incidence  approaches  90^, 
Also,  close  comparison  of  Figures  2A  and  2B  shows  a  slight 
flattening  of  the  distributions  as  the  slab  size  is  decreased^ 
all  of  which  points  to  the  fact  that  a  relatively  high  density 
of  collision  points  near  the  entrance  surface  results  in  a  flatten 
reflected  angular  distribution  than  a  more  uniform  distribution 
of  collision  points.  One  would  therefore  expect  that  a  decrease 
in  absorption  cross  section,  which  results  in  less  attenuation 
and  thus  a  more  uniform  collision  density,  would  result  in  more 
peaked  reflected  distributions.  This  effect  is  Illustrated 
in  Figures  3  and  4A,  B,  and  C  where  the  absorption  is  l/lOO  of 
the  total  cross  section.  In  this  case  the  reflected  distributions 
have  not  reached  equilibrium  even  at  10  mean  free  paths  (Fig.  3), 
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which  indicates  a  much  more  uniform  collision  density  in  the  slab 
than  was  the  case  with  a  higher  absorption  probability.  As  ex¬ 
pected;  Figures  4A  and  4B  show  the  reflected  angular  distributions 
to  be  more  peaked  toward  the  slab  normal  for  the  weak  absorber 
than  for  the  strong  absorber  (Figs.  2A  and  2B).  Figure  4C 
shows  the  effect  of  slab  size  on  the  reflected  distribution 
for  a  17.6°  Incident  angle.  This  effect  is  more  pronounced  for 
the  weak  absorber  than  for  the  strong  absorber  (not  shown).  As 
the  slab  size  is  Increased  the  reflected  distribution  increases 
over  the  entire  range  of  angles  but  more  rapidly  at  small  angles. 
The  result  is  that  the  angular  distribution  becomes  more  peaked 
as  the  slab  size  increases. 

The  effect  of  anisotropic  scatter  is  Illustrated  in  Figure  5. 
Because  the  scattering  cross  section  used  was  peaked  in  the  forward 
direction,  it  was  expected  that  the  collision  density  within  the 
slab  was  extremely  flat  for  near-normal  incidence,  thus  giving 
the  peaked  reflected  distributions.  For  grazing  incidence,  how¬ 
ever,  the  forward-peaked  cross  section  tends  to  concentrate  the 
collision  density  near  the  entrance  surface,  thus  giving  extremely 
flat  reflected  distributions.  In  fact,  in  the  extreme  case  of 
87.3°  incidence,  the  reflected  distribution  resembles  the  angular 
cross  section  because  much  of  the  reflection  is  due  to  singly 
scattered  particles. 
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2  X  ANGULAR  FLUX  (parficles/sec-unit  area-steradian) 


NPC  14,412 


PIGURi  4B.  REFLECnO  ANGULAR  FLUX  AS  A  FUNCTION  OF  REFLECTION 
ANGlEt  ABSORPTION  PROBABILITY  -  0.01 
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2  TT  X  ANGULAR  FLUX  (parficles/sec-unit  area-sferadian) 


NPC  14,413 


0  20  40  60  80  100  120  140 

REFLECTION  ANGLE  (deg) 

FIGURE  4C.  REFLECTED  ANGULAR  FLUX  AS  A  FUNCTION  OF  REFLECTION 
ANGLE:  ABSORPTION  PROBABILITY  ^  0.01 
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2  TT  X  ANGULAR  FLUX  (particles/sec-unit  area-s+eradian) 


NPC  14,414 


FIGURE  5.  REFLECTED  ANGULAR  FLUX  AS  A  FUNCTION  OF 
REFLECTION  ANGLE:  ANISOTROPIC  SCAHER 
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4,3.3  Transmission  Calculations 


For  the  transmission  calculations  shown  in  Figures  6  and  7 
it  was  assumed  that  scattering  was  isotropic  and  that  the  absorp¬ 
tion  cross  section  was  one-half  the  total  cross  section.  As  the 
slab  thickness  increases,  the  transmitted  distribution  for  thick 
slabs  decays  with  an  almost  constant  relaxation  length  -  as 
illustrated  in  Figure  6.  For  other  incident  angles  (not  shown), 
the  asymptotic  relaxation  length  is  the  same.  However,  cross 
plots  of  the  angular  flux,  such  as  Figure  7,  show  that  there  is 
still  some  tendency  for  the  angular  distribution  to  peak  in  the 
initial  direction  even  at  5.1  mean  free  paths.  This  is  probably 
because  the  high  absorption  cross  section  prevents  numerous 
scattering  collisions  which  would  tend  to  flatten  this  peaking. 
That  is,  the  particles  which  are  transmitted  have  suffered  very 
few,  if  any,  collisions. 

The  transmitted  angular  distributions  for  other  incident 
angles  are  similar  to  those  shown  in  Figure  7.  They  are  gen¬ 
erally  quite  peaked  toward  the  slab  normal.  This  is  explained 
by  the  fact  that  the  collision  density,  for  a  strong  absorber,  is 
much  greater  near  the  entrance  surface  and  thus,  as  mentioned 
earlier,  the  transmitted  distribution  is  determined  largely  by 
the  attenuation  suffered  along  the  slant  paths  from  the  high- 
collision— density  region  td  the  exit  point.  When  the  absorption 
cross  section  is  decreased  the  collision  density  inside  the  slab 
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is  5flatteuedj  resultiag  in  a  longer  asymptotic  relaxation  length 
(Fig.  8)  and  flatter  angular  distributions  (Fig.  9).  The  absorp¬ 
tion  cross  section  used  in  the  calculations  shown  in  Figures  8 
and  9  was  1%  of  the  total  cross  section.  The  introduction  of  a 
forward-peaked  scattering  cross  section  lengthens  the  relaxation 
length  still  more  (Fig.  10)  and  causes  a  very  slight  flattening 
of  the  transmitted  angular  flux  (Fig.  11), 

The  most  significant  result,  as  far  as  shielding  applications 
are  concerned,  is  the  asymptotic  behavior  of  the  transmission 
function.  In  all  cases  considered  here,  the  shape  of  the  trans¬ 
mitted  angular  distributions  approached  equilibrium  and  the 
magnitude  of  these  distributions  began  to  decay  exponentially 
with  a  relaxation  length  characteristic  of  the  material.  These 
results  and  their  obvious  applications  to  shielding  problems 
prompted  further  study  of  the  thick-slab  approximations  discussed 
in  Section  3.3.  In  these  approximations  the  transmission  function 
is  the  matrix  exponential 

Cx 

T(x)  -  e  , 

where  C  is  a  constant  matrix.  According  to  Sylvester's  Theorem, 

N  A.^x 

T(x)  -  E  e  z  , 
r-1  ^ 
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SLAB  THICKNESS  (mean  free  pafhs) 


FIGURE  6.  TRANSMiniD  ANGULAR  FLUX  AS  A  FUNCTION  OF  SLAB 
THICKNESS!  ABSORPTION  PROBABILITY  =  0.S 


2  TT  X  ANGULAR  FLUX  (particles/sec -unit  area- steradian) 


NPC  14,416 


0  20  40  60  80  100  120  140 


TRANSMISSION  ANGLE  (deg) 

FIGURE  7.  TRANSMIHED  ANGULAR  FLUX  AS  A  FUNCTION  OF  TRANSMISSION 
ANGLE:  ABSORPTION  PROBABILITY  ^  0.5 
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2  TT  X  ANGULAR  FLUX  (par+icles/sec-unitarea-steradian) 


NPC  14,417 


10-3  I _ I _ 1 _ I  ^ _ \ _ I 

0  2  4  6  B  10  12  14 

SLAB  THICKNESS  (mean  free  paths) 


FIGURE  8.  TRANSMiniD  ANGULAR  FLUX  AS  A  FUNCTION  OF  SLA8 
THICKNISSt  AiSORPTION  PROBAIILITY  =  0.01 
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X  ANGULAR  FLUX  (partlcles/sec-unit  area-sferadian) 


NPC  14,418 


FIGURE  9.  TRANSMiniD  ANGULAR  FLUX  AS  A  FUNCTION  OF  TRANSMISSION 
ANGLIt  ABSORPTION  PROBABILITY  =  0.01 
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2  TT  X  ANGULAR  FLUX  (parficles/sec- unit  area- s+eradian) 


SLAB  THICKNESS  (mean  free  paths) 


FIGURE  10.  TRANSMITTED  ANGULAR  FLUX  AS  A  FUNCTION  OF 
SLAR  THICKNESS:  ANISOTROPIC  SCAHER 


NPC  14,420 


TRANSMISSION  ANGLE  (deg) 

PIOURI  11.  TRANSMITTED  ANGULAR  FLUX  AS  A  FUNaiON  OF 
TRANSMISSION  ANOUt  ANISOTROMC  SCATHR 
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where  the  A.^,  are  the  N  eigenvalues  of  the  N  x  N  matrix  C  and  the 

are  the  matrices 
r 


2r  - 


g 

W  (Ag-X  ) 


r  /  s, 


where  U  is  the  Identity  matrix.  From  physical  considerations, 
it  is  seen  that  the  eigenvalues,  A^,  are  negative  (T(x)— ^0  as 
X— ^oo).  Also,  there  is  at  least  one  Am  such  that  |  Am  |— |  A^.! 
for  all  r  ^m  .  Thus,  for  x  large  enough,  the  mth  term  dominates 
and 

Am  X 

T(x)  e  zm  , 


which  is  of  the  form  indicated  by  the  numerical  results  of  this 
section.  That  is,  for  x  large  enough,  T(x)  decays  with  the 
relaxation  length  1/  |  Am  |  &nd  the  shape  of  the  angular  distribu<- 
tion  is  given  by  Zm  .  Thus,  from  a  knowledge  of  the  eigenvalues 
of  C,  one  can  determine  the  asymptotic  transmitted  angular  flux. 
The  accuracy  of  such  an  approximation  depends,  of  course,  on  the 
accuracy  of  the  thick— slab  approximation  and  the  equilibrium 
assumption  and  is,  therefore,  dependent  on  the  shield  material 
and  slab  thickness. 

The  fact  that  the  numerical  results  were  in  accord  with  the 
trends  predicted  from  collision-density  considerations  indicates 
that  this  approach  may  be  useful  in  further  studies  of  angular 
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distributions.  For  convenience,  the  conclusions  based  on  such 
considerations  and  verified  by  the  numerical  results  are  listed 
below. 


Perturbation 


Increase  absorption 
probabilities 

Increase  slab  thick¬ 
ness 

Increase  incident 
angle 

Increase  forward- 
scatter  probability 


Effect  on  Transmitted 
Angular  Distribution 


More  peaked 


More  peaked 


Little  change 
for  thick  slab 

Little  change 
for  thick  slab  — 
weak  absorber 


Effect  on  Reflected 
Angular  Distribution 

Flatter 


More  peaked 


Flatten^ 


More  peaked  for 
small  incident 
angle;  flatter 
for  grazing 
incidence 
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V.  CONCLUSIONS 


Although  these  initial  studies  indicate  that  definite 
advantages  may  be  gained  from  the  use  bf! 'imbedding,  techniques  in 
shielding  problems,  it  should  be  noted  that  regardless  of  the 
method  used  the  adequate  description  of  energy  and  angular  dis¬ 
tributions  requires  the  handling  of  large  amounts  of  data  and 
involves  a  considerable  number  of  computations.  Those  lamiliar 
with  the  limitations  of  present  computers  have  probably  wondered 
how  enough  Information  may  be  carried  in  the  machine  to  solve 
even  the  slab  problem  for  combined  neutron— gamma  transmission. 
This  may  not  be  possible  at  present.  Certainly  the  combination 
of  neutron,  gamma,  and  neutron-induced  gamma  penetration  calcula¬ 
tions  with  the  optimization  procedure  into  one  program  is  not 
feasible.  However,  the  author  believes  that  the  theoretical 
formulation  of  a  unified  approach  to  the  shielding  problem  is 
needed  -  regardless  of  its  immediate  practical  significance. 
Numerical  experiments  with  simplified  models,  such  as  the  one 
described  in  Section  IV,  may  uncover  new  problems  or  shed  some 
light  on  old  ones  -  they  may  even  result  in  the  discovery  of 
new  shielding  principles  which  are  obscured  by  our  present 
piece— meal  view  of  the  problem  and  our  attempts  to  be  realistic. 
Of  more  Immediate  importance,  perhaps,  is  the  possibility  that 
the  imbedding  formulation  of  the  shielding  problem  will  suggest 
new  approximations,  such  as  those  presented  in  Section  III. 
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Thus,  the  imbedding  method  offers  not  only  another  mathematical 
approach  to  conventional  shielding  problems  but  also  the  possibility 
of  the  development  of  new  physical  approximations  and,  ultimately, 
the  theoretical  formulation  of  the  shielding  problem  as  a  whole. 
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